Influence of correlated impurities on conductivity of graphene sheets: 
Time-dependent real-space Kubo approach 



T. M. Radchenko 
Solid State Electronics, Department of Science and Technology (ITN), 
Linkoping University, 60174 Norrkoping, Sweden and 
Deptartment of Solid State Theory, Institute for Metal Physics, 
NASU, 36 Acad. Vernadsky Blvd., 03680 Kyiv, Ukraine 

A. A. Shylau and I. V. Zozoulenko 

Solid State Electronics, Department of Science and Technology (ITN), 
Linkoping University, 60174 Norrkoping, Sweden 
(Dated: March 4, 2013) 

Exact numerical calculations of the conductivity of graphene sheets with random and correlated 
distributions of disorders have been performed using the time-dependent real-space Kubo formal- 
ism. The disorder was modeled by the long-range Gaussian potential describing screened charged 
impurities and by the short-range potential describing neutral adatoms both in the weak and strong 
scattering regime. Our central result is that correlation in the spatial distribution for the strong 
short-range scatterers and for the long-range Gaussian potential do not lead to any enhancement 
of the conductivity in comparison to the uncorrelated case. Our results strongly indicate that the 
temperature enhancement of the conductivity reported in the recent study (Yan and Fuhrer, Phys. 
Rev. Lett. 107, 206601 (2011)) and attributed to the effect of dopant correlations was most likely 
caused by other factors not related to the correlations in the scattering potential. 
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I. INTRODUCTION 

Investigation of charge transport in graphene and un- 
derstanding factors that affect its conductivity represent 
one of the central directions of graphene researchji~— This 
is motivated by both fundamental interest to graphene's 
transport properties as well as by potential applications 
of this novel material for electronics. It is commonly 
recognized that the major factors determining the elec- 
tron mobility in graphene are long-range charged impuri- 
ties trapped on the substrate and strong resonant short- 
range scatterers due to adatoms covalently bound to 
graphene i£ The nature of scatterers is directly reflected in 
the dependence of the conductivity on the electron den- 
sity, a = <r(n), and therefore investigation of this func- 
tion constitutes the focus of experimental and theoreti- 
cal research.^ For example, both the strong short-range 
potential^— and the long-range potential 7 -*^— lead to 
similar linear density dependence of the conductivity 
commonly observed in experiments Experiments 
also often show strong deviations from this linear depen- 
dence and an asymmetry with respect to the polarity 
of carriers, which can be attributed to a number of fac- 
tors such as the effect of realistic structural defects and 
impurities&i2r— , effect of contacts^—, effect of weak 
short-range impurities^^ - — , the ballistic or quasiballis- 
tic transport regime a 7 ' 27 " — and many others. 

Detailed studies of the density dependence of the con- 
ductivity in graphene often require exact numerical ap- 
proaches for transport calculations combined with ab 
initio calculations for microscopic properties of realistic 
scatterers Such the approaches can capture the es- 



sential features of the system at hand as well as can ad- 
dress transport regimes which are not accessible by other 
means. Exact numerical approaches can also be used to 
test a validity of conclusions of various semi-classical an- 
alytical approaches and applicability of approximations 
used in such the approaches. Among the most popu- 
lar numerical methods reported in the literature are the 
recursive Green's function technique^— and the time- 
dependent real space Kubo method££i2Ir— iiiil . The 
later method is especially suited to treat large graphene 
systems containing tens of millions of atoms with dimen- 
sions approaching realistic systems. 

Recently, Yan and Fuhrer— addressed the effect of cor- 
relation in the spatial distribution of disorder on the con- 
ductivity of graphene sheet by doping it with potassium 
atoms. They found that the conductivity of the system 
at hand increases as the temperature rises, and argued 
that this was caused by the enhancement of correlation 
between the potassium ions due to the Coulomb repul- 
sion. This conclusion, in turn, was based on the theoret- 
ical predictions of Li et al^ that the correlations in the 
position of long-range scatterers strongly enhances the 
conductivity. It should be noted that the approach of 
Li et alS^ is based on the standard Boltzmann approach 
within the Born approximation. However, the applica- 
bility of the Born approximation for graphene has been 
questioned in Rett where it was shown that predictions 
based on the standard semiclassical Boltzmann approach 
within the Born approximation for the case of the long- 
range Gaussian potential are in quantitative and quali- 
tative disagreement with the exact quantum-mechanical 
results in the parameter range corresponding to realistic 
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systems. (A discussion of various aspects of the applica- 
bility of the Born approximation for graphene and bilayer 
graphene can be found in Refsj£i2^£). It is therefore not 
clear to what extend the semiclassical predictions of Li 
et al££ based on the Born approximation are justified for 
the case of correlated impurities. Since the conclusions 
of the experiment of Yan and Fuhrer— rely essentially 
on the above semi-classical predictions, it is of interest 
to study the effect of spatial correlation of disorders by 
exact numerical methods. 

The main aim of the present paper is therefore to in- 
vestigate the influence of the spatial correlation of disor- 
der on the conductance of graphene sheets using an ex- 
act quantum mechanical approach. In the present study 
we utilize the time-dependent real-space quantum Kubo 
method£&^ ~ 23 ' 33 ~ — . allowing us to study large graphene 
sheets containing millions of atoms. We consider models 
of disorder appropriate for realistic impurities, including 
the Gaussian potential describing screened charged im- 
purities and the short-range potential describing neutral 
adatoms. The paper is organized as follows. The ba- 
sic model of the system under consideration including 
models for impurity potentials as well as the basics of 
the numerical Kubo method are formulated in Section 
II. Section III presents and discusses the obtained results 
and compares our numerical findings with available ex- 
perimental data and theoretical predictions. The conclu- 
sions are summarized in Section IV. A various details of 
the numerical approach are presented in the Appendix. 



II. TIGHT-BINDING MODEL AND 
TIME-DEPENDENT REAL-SPACE KUBO 
FORMALISM 
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FIG. 1. (Color online) A honeycomb graphene lattice. 

We model electron dynamics in graphene using the 
standard p-orbital nearest neighbor tight-binding Hamil- 



tonian defined on a honeycomb lattice 2 -! 3 -, Fig. [TJ 

h = -u^cta, +Y^v4a, (1) 

where cj and c, arc the standard creation and annihi- 
lation operators acting on a quasiparticle on the site 
i = (rn,n). The summation over i runs over the entire 
graphene lattice, while i' is restricted to the sites next to 
i; u = 2.7 eV is the hopping integral for the neighboring 
atoms i and i', and Vi is the on-site potential describing 
impurity scattering. 

In the present study we consider both short- and long- 
range impurities. The short-range impurities represent 
neutral adatoms on graphene and are modeled by the 
delta-function potential 

N imp 

Vi = E V i 5 H> ( 2 ) 

where Ni mp is the number of impurities on a graphene 
sheet. Estimations based on ab initio calculations and 
the T-matrix approach for common adatoms provide 
typical values for the on-site potential V 3 ■ — V) < 
80mi 9 i 19 ' 20 i 32 i 41 (For example, for hydrogen, Vo ~ 60u). 

The long-range impurities corresponds to charged ions 
situated typically on a surface of a dielectric substrate. 
We model them by the Gaussian potential commonly 
used in the literature, -^^^ 

^=E^exp(fc^), (3 ) 

where r, is the radius- vector of the site i, £ is the effective 
screening length, and the potential height is uniformly 
distributed in the range Uj G [—A, A] with A being the 
maximum potential height. 

We consider two cases of impurity distribution, 
random (uncorrelated) and correlated. In the first case 
of uncorrelated impurities the summation in Eqs. ([2"]l. 
is performed over the random distribution of impurities 
over the honeycomb lattice. In the second case impurities 
are no longer considered to be randomly distributed and 
to describe their spatial correlation we adopt a model 
used in Refi^ introducing the pair distribution function 
9(r), 

g(r) = { °' r < ro; (4) 
yv > \ 1, r > r , y ' 

where r is the distance between two impurities and the 
correlation length tq defines the minimum distance that 
can separate two impurities. (Note that for the ran- 
domly distributed (uncorrelated) impurities ro = 0). The 
largest distance ro max depends on the relative impurity 
concentration rif, the smaller the concentration n,, the 
larger ro max . Calculated values of ro max for the different 
relative impurity concentrations are presented in Table 
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FIG. 2. (Color online) A representative illustration of random 
and correlated distributions of impurities for the cases of (a), 
(b) short-range impurities, and (c), (d) Gaussian impurities 
for the impurity concentration rii = 2%. 
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TABLE I. Relation between the relative concentration of im- 
purities rii and the largest correlation distance ro max (ex- 
pressed in units of the carbon-carbon distance a — 0.142 nm). 



I, and representative examples of random and correlated 
distributions for the cases of the short- and long-range 
potentials are shown in Fig. 2. 

The transport properties of graphene sheets are cal- 
culated on the basis of the time-dependent real-space 
Kubo formalism where the dc conductivity a is extracted 
from the wave packet temporal dynamics governed by 
the time-dependent Schrodinger equation^— . This is a 
computationally efficient method scaling with a number 
of atoms in the system N, and thus allowing treating 
very large graphene sheets containing many millions of 
C atoms^^^i 

The calculation of the dc conductivity starts from the 
Kubo-Greenwood formula^ 



-Tr 



v x 8(E - H)v x 8{E - H) 



(5) 



where v x is the x-component of the velocity operator, E 
is the Fermi energy, VL is the area of the graphene sheet, 
and factor 2 accounts for the spin degeneracy. Let us 
introduce the mean quadratic spreading of a wave packet 

along the x-direction at the energy E, (AX 2 (E,t) 



X(t)-X(0)j ) where X(t) = W(t)XU{t) is the 

position operator in the Heisenberg representation, and 
U(t) = e ~ lHt / h is the time-evolution operator. The con- 
ductivity can then be expressed as the Einstein relation 

o = o xx = e 2 p(E) lim D(E,t), (6) 

t— >oo 

where p(E) = p/D, =Tr S(E - H) /fi is the den- 
sity of sates (DOS) per unit area (per spin), and the 
time-dependent diffusion coefficient D(E,t) is related to 
AX 2 (E,t)^ 



(AX 2 (E,t) 
D(E,t) = ± 



Tr 




f X H (t)-X(0)] 


2 ' 5{E-H) 


Tr 


6(E - H) 





(7a) 



(7b) 



Hence, calculation of a using the time-dependent real- 
space Kubo formalism requires numerical evaluation of 
the mean quadratic spreading of wave packets as pre- 
scribed by Eqs. © — (JZ]). Details of numerical calcula- 
tions of a and numerical solution of the time-dependent 
Schrodinger equation are given in the Appendix. 

The diffusion coefficient D(E,t) defined according to 
Eq. ([7]) is generally time-dependent and in different time 
intervals exhibits different temporal behavior depending 
on whether the transport regime is (i) ballistic, (ii) dif- 
fusive and (iii) localizedi 21 i 22 This is illustrated in Fig. [3] 
showing a temporal evolution of the diffusion coefficient 
for a graphene sheet with different impurity concentra- 
tions rii. The corresponding shapes of the wave packets 
for the different transport regimes are visualized in Fig. 
[3] for a some representative time t = 50 fs. In a system 
with no impurities (n^ = 0) electrons propagate ballisti- 
cally such that the mean spreading of the wave packet 

is J (^AX 2 (E, t)^ ~ vpt, where vp is the electron Fermi 

velocity. As a result, the diffusion coefficient increases 
linearly with time, D ~ vpt, with the slope being given 
by vp. In systems with a finite impurity concentration 
the ballistic regime lasts up to times Tbaii — di/vp where 
di = n i is the average distance between impurities. 

For times t > Tb a u the system enters the classical diffu- 
sive regime when the diffusion coefficient becomes inde- 
pendent of time. For the impurity concentration rii = 2% 
in Fig. |3] this regime starts at the time t w 25 fs when 
the diffusion coefficient D reaches its saturation. For 
larger times the diffusion coefficient starts to decrease 
due to quantum effects leading to the weak or the strong 
localization. This is in contrast to classical diffusion 
where D does not change when t — > oo. For example, 
for ni = 5% the localization effects becomes dominant 
already at t > 10 fs when the diffusion coefficient starts 
to decrease with time, see Fig. UJ (Note that for rii = 2% 
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FIG. 3. (Color online) (a) Temporal dependence of the dif- 
fusion coefficient for different concentration of strong short- 
range impurities. The on-site potential V ~ 37u; E = 0.2m. 

(b) — (e) Wave packet propagation in graphene lattice with 
different concentration of short-range impurities, (b) Initial 
random state of the size 1020x600 sites given by Eq. JOJ. 

(c) — (e) The wave packet shape after t — 50 fs for different 
impurity concentrations corresponding to ballistic, diffusive 
and localization regimes from (a). 

the decrease of D takes place at times exceeding the time 
interval shown in the figure). It should be stressed that 
in the present study we are interested in the diffusive 
transport regime when the diffusion coefficient reaches 
its maximum. Therefore, following Refsi 21 ' 22 , we replace 
in Eq. © lim^oo D{E,t) — > D max (E), such that the dc 
conductivity is defined as 

a = e 2 p(E)D max (E). (8) 

It is noteworthy that within the Boltzmann approach the 

2 

conductivity is given by UBoitz = e2 p(E) T ~i~i where r is 
the scattering time. Hence, it follows from Eq. (JSJ that 
the elastic length is related to the computed diffusion 
coefficient via l e = vpT — 2D max /vp. 

In most experiments the conductivity is measured as 
a function of the electron density n. Having calculated 
the DOS p(E) as described in the Appendix we compute 
the electron density n(E) = J_ p{E)dE — rti ons , where 

"ions = 3.9 x I0 15 cm~ 2 is the density of the positive ions 
in the graphene lattice compensating the negative charge 



of the p-electrons (note that for the ideal graphene lat- 
tice at the neutrality point n{E) = 0). (A dependence 
n = n(E) is illustrated in Fig. [4] for different scattering 
potentials). Combining the calculated n(E) with the con- 
ductivity cr(E) given by Eq. ©, we compute the density 
dependence of the conductivity, a = a(n). 

III. RESULTS AND DISCUSSION 

In this section we present results for the numerical con- 
ductivity of the graphene sheets with random and corre- 
lated impurities described by the short- and long-range 
potentials (Eqs. J2]) and (j3]) respectively). 

A. Short-range impurities 

Depending on the impurity strength, the conductivity 
of graphene sheets is expected to exhibit qualitatively dif- 
ferent density dependence for the weak and strong short- 
range scattering. The standard Boltzmann approach 
within the Born approximation predicts that the con- 
ductivity is independent on the electron density? 2 ^— ^ 



where Vf is the Fermi velocity in graphene, vfH = |ua, 
and Vo is the strength of the on-site potential in Eq. 
([2]). Vj = Vo- Going beyond the Born approximation 
one obtains that strong scatterers, with the accuracy 
up to logarithmic corrections, lead to a linear density 
dependency^ - — 

v=^r— (lnV^o) , (10) 

where i?o is the scatterer's radius. Note that the condi- 
tion for the validity of the Born approximation in Eq. ^ , 
can be presented in the form^^ (|Vb|/u)~ <C y^lnrie, 
where n e is the relative electron density (i.e. the elec- 
tron density per C atom). For realistic electron densi- 
ties, n e < 0.01, this condition loosely corresponds to |Vb| 
< u, which we adopt as a definition of the weak impu- 
rity scattering. The condition \Vq\ S> u (when Eq. (JTUJ) 
is expected to be valid), defines the case of the strong 
impurity scattering. Note that most of adatoms present 
in graphene falls into the second case of the strong scat- 
tering with the effective on-site potentials in the range 
|Vb|/u« 10-80^1 

1. Short-range potential, strong scattering, \Va\ S> u 

Figures @] (a), (d) show the DOS and the electron den- 
sity n = n(E) in a graphene sheet with strong scatterers, 
Vb = 37m and the concentration rij = 2%. The calcu- 
lated dependencies are practically identical for the cases 
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FIG. 4. (Color online) Density of states (DOS) and the relative charge carrier concentration n e (the number of electrons per 
C-atom) as function of the energy E for (a) short-range strong, (b) short-range weak, and (c) Gaussian potential for random 
and correlated impurities, (c) Dependencies n = n(E) for short- and long-range potentials. 



of random and correlated distributions of impurities. The 
calculated DOS shows a pronounced peak in the vicinity 
of the Dirac point, —0.3 < — < 0.1 which reflects a for- 
mation of the impurity band as discussed in Ref£. Note 
that the peak is shifted with respect to E = which 
is related to the asymmetry of the impurity potential 
(Vq > 0). The shift of the impurity band also leads to 
the asymmetry in n = n(E) dependence, see Fig. 0](d). 

Figure [5] (a) shows the time-dependent diffusion coeffi- 
cient for different impurity concentrations for the case of 
random distribution of impurities. The temporal depen- 
dencies D = D(t) demonstrate that the system at hand 
reaches and stays in the diffusive transport regime for the 
impurity concentration n.j = 1% and 2% at times t > 70 
and 30 fs correspondingly. For the impurity concentra- 
tion m = 0.5% the diffusive regime is reached at times 
outside those shown in the figure (t > 130 fs), whereas 
for n, = 5% the system almost immediately reaches the 
localization regime. 

To analyze the density dependence of the conductivity 
we choose a representative concentration rij = 2% and 
t = 80 fs when D{t) exhibits a saturation corresponding 
to a well-defined diffusive regime, see Fig. 0(b). The ran- 



domly distributed impurities show a quasi-linear density 
dependence of the conductivity. It is worth to note that 
the calculated conductivity fully reproduces numerical 
results reported by Yuan et al& who used a similar time- 
dependent Kubo approach. The corresponding theoret- 
ical prediction exhibits pronounced deviations from the 
linear dependence caused by logarithmic corrections in 
Eq. (fTT)]) . These deviations are much stronger than those 
typically seen experimentally^ This is related to the fact 
that in realistic experimental samples the electron densi- 
ties are lower than those used in our calculations^ such 
that at smaller densities Eq. ([TU]) also exhibits a quasi- 
linear dependence. In the vicinity of the Dirac point the 
conductivity flattens out which is related to a transport 
regime due to a formation of the impurity band. Ap- 
parently, the theoretical prediction Eq. (TTU|) docs not re- 
produce this transport regime. For larger densities away 
from the Dirac point the calculated Kubo conductivity 
differs by a factor of ~ 2 from the theoretical predictions 
given by Eq. ([TU)) . Note that the calculated dependence 
a = cr(n) is shifted with respect ton = 0, which is caused 
by the asymmetry in the dependence n = n(E) due to 
the impurity band as discussed above. 

Figure O (b) also shows the dependence a = <x(n) 
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short-range strong 




FIG. 5. (Color online) (a) Diffusion coefficient D = D(t) for different concentration of the short-range strong impurities; 
Vo = 37u, E — 0.2u. (b) The conductivity a as a function of the relative electron density n e (the number of electrons per 
C-atom) for random and correlated impurities, m = 2%. The dot-dashed line corresponds to Eq. (10), which is shifted to the 
charge neutrality point at n e ~ 0.02. 



for the case of correlated impurities with the correlation 
lengths ro = 0.5ro m ax and r$ = ro max .The central result 
is that correlation in the impurity distribution practically 
does not affect the conductivity of the system even for 
the largest correlation length rg max • 



2. Short-range potential, weak scattering, Vo <u- 

Let us now turn to the case of weak scattering. We 
consider two cases, the symmetric potential, Vq = ±it, 
and the asymmetric one, Vo = u. Figures 2] (b) , (d) show 
the DOS and the electron density n — n(E) for the case 
of rii = 2%. As in the case of the strong scatterers, the 
calculated dependencies are practically identical for ran- 
dom and correlated distributions of impurities. However, 
in contrast to the case of the strong scatterers, the DOS 
does not form the impurity band close to the Dirac point. 

Figure [6] shows the time-dependent diffusion coeffi- 
cients D(t) and the density dependence of the conduc- 
tivity for rii = 3% with random and correlated distribu- 
tions of disorders. In contrast to the case of the strong 
scattering potential, in the present case the diffusion co- 
efficient reaches its maximum at the same time t m 120 fs 
independent of the impurity concentration. The density 
dependence of the conductivity, a = cr(n), for symmet- 
ric scatterers is shown in Fig. [6] (b). In agreement with 
the Boltzmann predictions, Eq. (O, the calculated Kubo 
conductivity is independent of the carrier density, even 
though its calculated value differs by a factor of ~ 3 from 
the corresponding theoretical predictions. The calculated 
a is independent on n for all densities except those in the 
vicinity of the Dirac point where the conductivity shows 
a pronounced dip. This can be explained by the fact 
that close to the Dirac point the graphene sheets is in the 
electron- hole puddle density regime with strong potential 



variations comparable to the Fermi energy of electrons. 
Because Eq. © is not valid in such the regime, the cal- 
culated density strongly deviates from the predictions of 
Eq. (JSJ • Figure \§\ (b) also shows the the conductivity for 
the case of correlated impurities with ro = 0.6ro m ax and 
ro = romax- Similarly to the case of strong scatterers the 
correlation in the impurity distribution practically does 
not affect the conductivity of the system at hand. 

The density dependence the conductivity, a = cr(n), 
for the case of asymmetric potential, Vo = u, is presented 
in Fig. [S] (d). It has two features that are different from 
the case of the symmetric potential. First, for negative 
electron densities, the calculated conductivity is indepen- 
dent of n, which is consistent with the Boltzmann pre- 
dictions, Eq. ©. However, for positive densities, the 
conductivity shows a sublinear density dependence dis- 
tinct from Eq. ([9]). Second, for the case of correlated 
impurities the conductivity is increased by up to 30% in 
the region of negative densities (when a wconst), whereas 
a is not affected by the correlations for positive densities 
where the behavior is sublinear. It should be stressed 
that the Boltzmann theory predicts the same density de- 
pendence regardless whether the potential Vo symmetric 
or not, whereas our numerical calculations show a clear 
difference between these two cases. Note that calcula- 
tions with the asymmetric potential, Vo = —it, show the 
density dependence of the conductivity which is mirror- 
symmetric to the case of Vq = u. 



B. Long-range Gaussian impurities 



Let us now turn to the case of the long-range Gaussian 
potential, Eq. ([3]). The graphene conductivity calculated 
within the the standard Boltzmann approach in the Born 
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FIG. 6. (Color online) (a), (c) The diffusion coefficient D — D(t) for for short-range weak impurities with different concentrations 
rii', E — 0.2m. (b), (d) The conductivity a as functions of the relative electron density n e (the number of electrons per C-atom) 
for random and correlated short-range weak impurities, rii — 3%. Dot-dashed lines shows the Boltzmann predictions according 
to Eq. ©. Panels (a), (b) correspond to the symmetric potential, Vo = ±u, and (b), (d) for asymmetric one, Vo = u. 



approximation reads ^ 42 ' 46 , 

_ 4e 2 nn^e 7 " 1 ^ f const, \z\ < 1, 
°~ h Kh{ime) ^ \ n 3 / 2 , \z\ » 1. {) 

where z = mi£ 2 = (^j^j, h is the modified Bessel 

function and K w 40.5n ?; (A/it) 2 (^/V3a) 4 with a be- 
ing the carbon-carbon distance and A being the Fermi 
wavelength. The condition \z\ <C 1 describe the case of 
quantum scattering when the Fermi wavelength is larger 
than the effective screening length, A ^> £, while the 
opposite condition \z\ ^> 1 corresponds to the case of 
classical scattering, A « ^. Because the semiclassical 
approach predicts two distinct regimes of conductivity, 
in our numerical calculations we explore the parameter 
space corresponding to these two regimes. We consider 
£ = 4a and £ = 16a, which in the considered density in- 
terval | n e | < 0.06 correspond respectively to \z\ < 2, and 
\z\ < 35. Note that the regime \z\ < 1 is appropriate to re- 



alistic experimental samples, whereas the carrier density 
\z\ > 1 is too high to be achieved in the experimental 4 ^ 

Figures 0] (c), (d) show the DOS and the electron den- 
sity n = n(E) for the case of rii = 2%. The Gaussian 
potential significantly smooths the density of states es- 
pecially in the region of the van-Hove singularities. As in 
the case of the short-range scatterers, the calculated de- 
pendencies are practically identical for random and cor- 
related distributions of impurities. 

Figure [7| shows the time-dependent diffusion coeffi- 
cients D(t) and the conductivity a = cr(n) for n, = 4% 
with random and correlated distributions of disorders. 
As in the case of the weak short-range disorder in the 
present case the diffusion coefficient reaches its maximum 
value D max at the same time for all impurity concentra- 
tions studied 1% < rii < 5% (t = 130 fs respectively 110 
fs for £ = 4a respectively 16a). 

The density dependence of the conductivity for £ = 4a 
is shown in Fig. [7] (b) in the density interval \z\ < 2. 
The conductivity exhibits the linear density dependence 
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FIG. 7. (Color online) (a), (c) The diffusion coefficient D = D(t) for the Gaussian impurities with different concentrations rn\ 
E — 0.2u. (b), (d) The conductivity a as functions of the relative electron density n e (the number of electrons per C-atom) 
for random and correlated Gaussian impurities, m — 4%. A dot-dashed line shows the Boltzmann predictions according to Eq. 
(|lip . Inset in (d) is plotted in the logarithmic scale. Panels (a), (b) and (c), (d) correspond to £ = 4a and £ = 16a respectively. 



which well extends into the classical regime \z\ > 1. 
This is in agreement with previous numerical calculations 
performed using the Green's function technique^^S, also 
showing linear (or quasi-linear) density dependence of the 
conductivity. On the contrary, the standard Boltzmann 
approach predicts that for \z\ <C 1 the conductivity is 
independent on the electron density, a =const. For the 
present case of the Gaussian potential it has been ar- 
gued that the standard Boltzmann approach (leading to 
a =const) is not valid because the Born approximation 
is not well justified^ 

Figure [Jj (d) shows the conductivity for the case of 
Gaussian disorder for £ = 16a in the density interval 
\z\ < 35. Even in this case the calculated conductivity 
shows the linear density dependence, a cx n. For \z\ ^> 1 
the Boltzmann approach predicts the superlinear depen- 
dence a oc rt 3 / 2 (this is clearly discernible in the loga- 
rithmic scale as shown in the inset of Fig. [JJ (d)). We 
however do not reach such high densities to reproduce 
this asymptote. (This regime of the high densities was 
analyzed in Ref^ 2 -). 



Figures [7J (b) and (d) also shows the conductivity for 
the case of the correlated impurity distribution for £ = 4a 
and £ = 16a. Even for the maximal correlation between 
impurities, tq — ro ma x the correlation practically does 
not affect the conductivity. This represents the central 
result of the present section. 

C. Comparison to the experiment 

Let us now compare our numerical data with avail- 
able experimental results and theoretical predictions. Re- 
cently, Yan and Fuhrcr reported conductivity measure- 
ments of potassium doped graphene^. They found that 
with the increase of the temperature the conductivity of 
doped graphenc increases by up to a factor of 2. They 
attributed this effect to the enhancement of the spatial 
correlation between potassium ions due to the mutual 
Coulomb repulsion which, according to the recent the- 
ory of Li et aL— leads to the conductivity enhancement. 
Our exact numerical calculations however do not sup- 
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port this conclusion. Indeed, we demonstrated that the 
spatial correlations of charged impurities (modeled by 
the long-range Gaussian potential) has no effect on the 
conductivity of graphene sheets. Our results therefore 
suggest that the enhancement of the conductivity with 
increase of the temperature in experiments of Yan and 
Fuhrer— is most likely caused by other factors not related 
to the correlations of impurities. Moreover, as no direct 
evidence of the spatial correlation of the potassium ions 
was presented in the above article, it is not clear whether 
impurities in experimental samples are correlated at the 
first place, and whether this correlation (if any) increases 
with the temperature. 

Our exact numerical results do not apparently support 
theoretical predictions of Li et al£2- that the correlations 
in the impurity positions lead to the enhancement of the 
conductivity. As mentioned before, it has been shown 
that for the long-range Gaussian potential in a parame- 
ter range corresponding to realistic systems the standard 
Boltzmann predictions are in quantitative and qualita- 
tive disagreement with the exact numerical results^ The 
reason for that is the utilization of the Born approxi- 
mation which relies on the unperturbed wave functions 
of ideal graphene without impurities. This is justified 
only for the case of weak scattering and apparently fails 
for the long-range Gaussian scatterers in the parameter 
range corresponding to realistic systems. (For the discus- 
sion of the Born approximation for graphene and bilayer 
graphene see Refsji^^). This explains the discrepancy 
the between the exact results and predictions of Li et al£2- 
which are essentially based on the standard Boltzmann 
approach in the Born approximation. 



only when the potential is asymmetric, i.e. V — Vq or 
V = —Vq. No enhancement of the conductivity is found 
for the symmetric weak short-range potential, V = ±Vq. 

Using our results we analyze the recent experiment of 
Yan and Fuhre r 3839 where the temperature increase of 
the conductivity was attributed to the enhancement in 
the spatial correlation of the adsorbed potassium ions. 
Our numerical findings do not sustain this interpretation 
and our results strongly suggest that the enhancement 
of the conductivity reported in the above study is most 
likely caused by other factors not related to the correla- 
tions of impurities. 

Our numerical calculations do not support theoreti- 
cal predictions of Li et al££ that the correlations in the 
impurity positions for the long-range potential lead to 
the enhancement of the conductivity. We attribute this 
to the utilization of the standard Boltzmann approach 
within the Born approximation which is not justified for 
the case of a long-range potential in the parameter range 
corresponding to realistic systems. 
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Appendix: Numerical calculation of the dc 
conductivity a on the basis of the time-dependent 
real space Kubo method 



IV. CONCLUSIONS 

Using an efficient time-dependent real space Kubo for- 
malism we performed numerical studies of conductivity 
of large graphene sheets with random and correlated 
distribution of disorder. In order to describe realistic 
disorder we used models of the short-range scattering 
potential (appropriate for adatoms covalently bound to 
graphene) and the long-range Gaussian potential (ap- 
propriate for screened charged impurities on graphene 
and/or dielectric surface). The calculations for the uncor- 
relatcd potentials are compared to the corresponding pre- 
dictions based on the semiclassical Boltzmann approach 
and to exact numerical calculations performed by differ- 
ent methods. 

We find that for the most important, experimentally 
relevant cases of disorder, namely the strong short-range 
potential and the long-range Gaussian potential, the cor- 
relation in the distribution of disorder docs not affect the 
conductivity of the graphene sheets as compared to the 
case when disorder is distributed randomly. This rep- 
resent the main result of our study. We find that the 
correlations lead to the enhancement of the conductivity 
only for the case of the weak short-range potential and 



1. Calculation of the DOS p(E) and the 
time-dependent diffusion coefficient D(E, t) 

Numerical calculation of the dc conductivity a requires 



=Tr 



S(E - H) 



and the 
Let 



computation of the DOS p(E) 

time-dependent diffusion coefficient D(E, t), Eq. ([7J 
us start with an algorithm for calculation of p(E). 

Express the density of states p{E) as a sum over the 
local densities of states p(E) = Pi{E), where the 
summation is performed over all the sites of graphene 
lattice N. In its turn, the local density of states is given 
by the imaginary part of the diagonal elements of the 
Greens function, pi(E) = —-^sGii(E + ie), where a small 
e — > is introduced in order to smooth the peaks in the 
DOS. 

There is an efficient numerical algorithm for calcula- 
tion of the diagonal elements Ga . It first starts with the 
tridiagonalization procedure when the Hamiltonian is re- 
duced to the tridiagonalized form by passing to a new 
basis. Then the first diagonal element of the Green's 
function, Gu, is computed using the standard continued 
fraction technique. The details of this algorithm are pre- 
sented in Appendix B. Note that this algorithm scales 
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as the number of sites N, with the most time-consuming 
part being the triagonalization procedure. Having calcu- 
lated the local density of states for the first site, pi(E), 
one can in principle repeat similar calculations for all re- 
maining N — 1 sites in order to recover the total DOS 
p{E). This would apparently require total computation 
efforts that scale as ./V 2 . For a disordered systems treated 
in the present paper one can adopt a different way to eval- 
uate p{E) that still keeps the total computation efforts 
of the order of N. It relies on the fact that a sufficiently 
large subsystem of the total system has the same density 
of states as the original one. Our computation of p(E) 
proceeds then as follows. 

Choose a subsystem occupying M sites of the graphene 
lattice. Construct a random state \ip r ) extended over 
these M sites, 



e 2i7Tai \i) 



(A.l) 



with ati being the random phase in the interval [0,1], 
\i) = c\ |0), and the summation is extended over the sites 
of the chosen subsystem. An example of such a random 
state is shown in Fig. [3](b). Transform the original tight- 
binding Hamiltonian Eq. (TTJ defined in the basis {r{\ by 
passing on to a new basis as follows. Choose the first ele- 
ment of the new basis as |1} = iV^-an), where |^ r an) is the 
random state defined according to Eq. 



(TA~Tj) . Tridiag- 
onalize the Hamiltonian as prescribed in Appendix and 
calculate pi{E) = —^5sGu(E + is) using the continued 
fraction technique. Repeat this procedure, if needed, for 
different distributions of disorder and average p\ (E) over 
these disorder realizations. The calculated value of p\{E) 
corresponds to the DOS per carbon atom of the system 
at hand. Note that calculation of the remaining N — 1 
matrix elements Gu is not needed such that the total 
computational efforts still scale as N. 

Calculation of the trace 



Tr 



(x H {t)-x(o) 



S(E - H) 



in the expression 
performed in a similar 



for D(E,t) in Eq. 
way. It can be shown^ that this trace can be rep- 
resented as a sum of the local densities of states 



Tr 



(x H (t) - X(of) S(E 
functions $i(t), where 



H) 



Ef Pi(E,t) for the 



N 



i=l 



p i (E,t) = Y^(MtME-H)\%(t) 



(A.2) 



and 



|*i(t)} = XU\^i) - U (l|Vi>) , (A.3) 

where X is the position operator (i.e. the ^-coordinate), 
U is the time evolution operator and {ipi} is the orthog- 
onal basis set (corresponding to e.g. the eigenstates of 
the Hamiltonian Eq. (fT])). Next, we pass on to the new 



basis setting its first element |1} = ^(f)), where \^>i(t)) 
is given by Eq. (|A.3[) where \ipi) = |V>ran) is chosen as 
the random state defined by Eq. (jA.ll) . Time evolution 
of the random initial state in Eq. (|A.3[) is calculated by 
means of the Chebyshev method as described below. We 
then reduce the Hamiltonian to the tridiagonalized form 
and use the continued fraction technique to calculate lo- 
cal density of states pi(E, t) in the new bases. As argued 
above, for a sufficiently large random state, this local 
density of state well approximates the density of states 
of the whole system. 



2. Exact solution of the time-dependent 
Schrodinger equation by means of the Chebyshev 
method 



In this appendix we present an efficient method for so- 
lution of the time-dependent Schrodinger equation based 
on the expansion of the time evolution operator in an 
orthogonal set of Chebyshev polynomials. 

We start from the time-dependent Schrodinger equa- 
tion complemented by the initial condition 



HW)) 



m = o)) = \,M 



(A.4) 



with H being the tight-binding time-independent Hamil- 
tonian, Eq. (fT]). Its formal solution can be expressed via 
the time-evolution operator U(t), 



(A.5) 



In order to expand U (t) in a set of the Chebyshev poly- 
nomials T n {x) (which are defined in the interval x € 
[— 1 ; 1] ) , we we first renormalize the Hamiltonian such 
that its spectrum lies in the above interval, 



2H - (E n 



E m i n )I 



EL 



E„ 



(A.6) 



where -E max and -E m i n are the largest and the smallest 
eigenvalues of the original Hamiltonian Eq. (TTJ. (In or- 
der to calculate -E max and ^min we use a computational 
routine that estimates the largest/smallest eigenvalues of 
the operator without calculation of all the eigenvalues). 

Expanding U(t) in Chebyshev polynomials in Eq. 
(|A.5p we obtain for the wave function, 



(A.7) 



n=0 



where the functions |$ n ) = T n (H noTm ) |i/>o) are calcu- 
lated using the recurrence relations for the Chebyshev 
polynomials, 



1$, 



l*n-l) 



(A. 
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with |$q) = To(H noim ) \ipo) — \ipo)- The expansion co- 
efficients c n (t) are calculated making use of the orthogo- 
nality relation for the Chcbyshev polynomials 



c n (t) = 2e~ 



t +g min )t 



{-i) n Jn 



^max E n 

2h 



-t 



(A.9) 

For large t the expansion coefficients c n (t) become expo- 
nentially small. This leads to the fast convergence of the 
expansion series Eq. (|A.7[) . an d makes the Chebyshev 
method very efficient for calculation of the temporal dy- 
namics. For instance, in our simulations we take 10000 
iterations for 300 fs of elapsed time. 



3. Continued fraction technique and 
tridiagonalization of the Hamiltonian matrix 

a. Continued fraction technique 



Consider a Hamiltonian matrix in a tridiagonal form, 

\ 



/ ai j3i 



tri 



Pi 





a 2 



V 



o 



«3 Pa 



Av-i 

Ot N J 



(A.10) 



The continued fraction technique provides an efficient 
way to calculate the first diagonal element Gxi of the 
Greens function G = (EI — iftri) -1 without a need for 
computing the whole Greens function and all the eigen- 
values/eigenfunctions of the Hamiltonian. 

Let us denote A; = Gu. From the definition of the 
Greens function we obtain, 

(E-afjXi-Pi-xXi-x-foXi+x = 0, 2 < i < N-l (A.ll) 

with (E — aN)^N — Pn-iXn-i = and (E — Oi)Ai — 
PxX% = 1. Expressing sequentially Xn by Ajv_i, Aat_i via 
Aat_2 and Ajv-3 we express Gxx as a continued fraction: 



1 



E — on 



(A.12) 



E-OL2-- 



In the above equation we truncated the order- N contin- 
ued fraction at the fraction M < N by introducing the 
self-energy T*(E), 



£(£) = 



1 



1 



E - a M 



PI 



E-a M - P 2 M E{E) 



E — O.M- 



that includes all the remaining terms M + 1 < 
Solving Eq. (|A.13|) one easily obtains 



£(£) 



E-O.M- W0m ~{E- a M ) 



2(3 2 



< N. 



(A.14) 



M 



The number of terms M included in the summation in 
Eq. (|A.12|) is determined from the condition for the con- 
vergence of Gu. For instance, in our calculations with 
the N x N Hamiltonian matrixes with = 1.7 x 10 6 and 
6.8 x 10 6 it is sufficient to choose M m 2000 and 4000 
respectively. 



b. Tridiagonalization of the Hamiltonian matrix 

In order to utilize the continued fraction technique to 
calculate Gu the Hamiltonian should be transformed to 
the tridiagonalizcd form, Eq. (|A.10|) . This is done by 
constructing a new orthogonal basis as described below. 
We start by selecting the first basis vector |1}. We use 
curled brackets |i} to denote the new basis vectors and 
straight brackets \i) to denote the old ones. If the tridi- 
agolization is performed in order to find the local density 
if states on the i-th site of the system at hand, the first 
basis vector is selected as |1} = \i), where |i) = c-|0). In 
most our calculations we select the first basis vector as 
the random state occupying M sites, |1} = |'0ran)j with 
|Vw) being given by Eq. (|A.1|) . 

We require that the Hamiltonian in the new basis be of 
the tridiagonal form (|A.10[) . By operating H\i} we arrive 
to the following equations, 

fl-|l} = ai|l}+0i|2}, (A.15a) 

H\i} = ft_i|i - 1} + cti\i} + /3i\i + 1}, 2 < i < (SL.±Sb) 

H\N}^/3 N ^\N-l} + a N \N}. (A.15c) 

Using Eq. (|A.15ap and the orthogonality relation {1|2} = 
we obtain the second basis vector and the matrix ele- 
ments oti and /3i, 

\2} = -±=(H\l}-a 1 \l}), (A.16) 



/C2 

ai={l|^|l}, px={2\H\l], 

where the normalization coefficient G2 (as well as all 
other normalization coefficients G,,2 < i < N) are ob- 
tained from the normalization requirement {i\i} = 1. 

We then proceed to Eq. (|A.15b|) and recursively calcu- 
late the basis vectors \i},2 < i < A — 2, and correspond- 
ing matrix elements ai and Pi, 

\i + 1} = — ^= (h\%} - pi-x\i - 1} - Oil*}) , (A.17) 

a t = {i\H\i}, pi +1 ={1 + 1\H\1], 2 < 1 < N - 2 
Finally, from Eq. (|A.15c|) we obtain 
1 



\N} 



(H\N} - a N \N}^ , a N = {N\H\N}, 



(A. 13) which concludes the tridiagonalization procedure. 
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FIG. 8. Propagation of wave packets of two different initial sizes: 170x100 and 1020x600 on a graphene sheet with random 
distribution of the short range strong impurities V = 37u, m = 1%. 



4. Role of the initial wave packet size in the 
averaging over impurity configurations 

For a given configuration of impurities, the conductiv- 
ity of the system at hand can be sensitive to details of 
the potential and thus can vary from one impurity real- 
ization to another. The conductivity therefore should be 
averaged by e.g. performing many calculations for dif- 
ferent impurity distributions. Within the present time- 
dependent Kubo approach the averaging can be done in 
a much more efficient way by simply increasing a size of 
the initial random state without need for many different 
calculations for different impurity realizations. In this 
section we investigate how one can achieve an efficient 
averaging of the conductivity by varying the size of the 
wave packet. 

Figure [5] shows a temporal evolution of two random ini- 
tial states l^ran) (Eq. (|A.1|) ) of different sizes, 170x100 
and 1020x600 respectively. The corresponding time- 
dependent diffusion coefficients D(t) are shown in Fig. 
[5] for different electron energies E. For the case of 
the smaller packet the calculated temporal dependen- 
cies show fluctuations caused by the interference within 
the wave packet. In contrast, for the case of the larger 
packet these fluctuations are efficiently averaged out and 
D(t) exhibit a smooth monotonic behavior for all con- 



sidered energies. These self-averaging features of the 
wave packets temporal dynamics manifest themselves in 
the density dependencies of the conductivity a = a{n). 
Figure [TU] (a),(b) show the dependencies a = <j(n) for 
the above wave packets for different impurity realizations 
(Note that the concentration of impurities in all cases 
is the same.) For the case of the smaller wave packet the 
conductivity exhibits significant fluctuations, whereas for 
the case of the larger wave packet the conductivity curves 
practically coincide. It is important to stress that the av- 
eraged values of the conductivity are the same in both 
cases, see Fig. [TU] (c). Thus, in most our calculations 
we choose the wave packet to be sufficiently large (typ- 
ically 1020x600) such that no additional averaging over 
impurity realization is needed. 

It is noteworthy that a larger wave packet apparently 
reaches the boundary of the computational domain ear- 
lier than the smaller one, c.f. Fig. [8] (c) and (f). The 
size of the computational domain should be therefore 
sufficiently large, such that the maximum value of the 
diffusion coefficient, Eq. (|7a[) . is reached before the 
wave packet hits the boundaries. In our calculations we 
used the honeycomb lattices of the sizes 1700x1000 and 
3400 x 2000 sites (corresponding to 2 10 x 2 1 and 420 x 420 
nm 2 respectively). 
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FIG. 9. (Color online) Time-dependent diffusion coefficient D(t) at different energies E for the wave packets (WP) of Fig. [8] 




FIG. 10. (Color online) (a), (b) Conductivity vs. energy for the wave packet (WP) of Fig. [5] for ten different impurity 
configurations, (c) Averaged values of conductivities for the wave packets from (a) and (b). 
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